Position - Linearization¶

Overview¶

Developer Note: if you may make a PR in the future, be sure to copy this notebook, and use the gitignore prefix temp to avoid future conflicts.

This is one notebook in a multi-part series on Spyglass.

To set up your Spyglass environment and database, see the Setup notebook
For additional info on DataJoint syntax, including table definitions and inserts, see the Insert Data notebook

This pipeline takes 2D position data from the IntervalPositionInfo table and "linearizes" it to 1D position. If you haven't already done so, please generate input data with either the Trodes or DLC notebooks (1, 2, 3).

Imports¶

In [ ]:

Copied!

%reload_ext autoreload
%autoreload 2
%reload_ext autoreload
%autoreload 2

In [ ]:

Copied!





import os
import pynwb
import numpy as np
import datajoint as dj
import matplotlib.pyplot as plt

# change to the upper level folder to detect dj_local_conf.json
if os.path.basename(os.getcwd()) == "notebooks":
    os.chdir("..")
dj.config.load("dj_local_conf.json")  # load config for database connection info

import spyglass.common as sgc
import spyglass.position.v1 as sgp

# ignore datajoint+jupyter async warnings
import warnings

warnings.simplefilter("ignore", category=DeprecationWarning)
warnings.simplefilter("ignore", category=ResourceWarning)
import os
import pynwb
import numpy as np
import datajoint as dj
import matplotlib.pyplot as plt

# change to the upper level folder to detect dj_local_conf.json
if os.path.basename(os.getcwd()) == "notebooks":
    os.chdir("..")
dj.config.load("dj_local_conf.json")  # load config for database connection info

import spyglass.common as sgc
import spyglass.position.v1 as sgp

# ignore datajoint+jupyter async warnings
import warnings

warnings.simplefilter("ignore", category=DeprecationWarning)
warnings.simplefilter("ignore", category=ResourceWarning)

[2023-07-28 14:45:50,776][INFO]: Connecting root@localhost:3306
[2023-07-28 14:45:50,804][INFO]: Connected root@localhost:3306

Retrieve 2D position¶

To retrieve 2D position data, we'll specify an nwb file, a position time interval, and the set of parameters used to compute the position info.

In [2]:

Copied!

nwb_file_name = "chimi20200216_new.nwb"
nwb_copy_file_name = sgc.nwb_helper_fn.get_nwb_copy_filename(nwb_file_name)
nwb_copy_file_name
nwb_file_name = "chimi20200216_new.nwb"
nwb_copy_file_name = sgc.nwb_helper_fn.get_nwb_copy_filename(nwb_file_name)
nwb_copy_file_name

[2022-08-04 16:16:35,902][INFO]: Connecting zoldello@lmf-db.cin.ucsf.edu:3306
[2022-08-04 16:16:35,953][INFO]: Connected zoldello@lmf-db.cin.ucsf.edu:3306

/home/zoldello/anaconda3/envs/spyglass/lib/python3.9/site-packages/position_tools/core.py:3: DeprecationWarning: Please use `gaussian_filter1d` from the `scipy.ndimage` namespace, the `scipy.ndimage.filters` namespace is deprecated.
  from scipy.ndimage.filters import gaussian_filter1d

Out[2]:

'chimi20200216_new_.nwb'

We will fetch the pandas dataframe from the IntervalPositionInfo.

In [3]:

Copied!





position_info = (
    sgc.common_position.IntervalPositionInfo()
    & {
        "nwb_file_name": nwb_copy_file_name,
        "interval_list_name": "pos 1 valid times",
        "position_info_param_name": "default",
    }
).fetch1_dataframe()
position_info
position_info = (
    sgc.common_position.IntervalPositionInfo()
    & {
        "nwb_file_name": nwb_copy_file_name,
        "interval_list_name": "pos 1 valid times",
        "position_info_param_name": "default",
    }
).fetch1_dataframe()
position_info

/stelmo/nwb/analysis/chimi20200216_new_ZF24K2JHC1.nwb

Out[3]:

	head_position_x	head_position_y	head_orientation	head_velocity_x	head_velocity_y	head_speed
time
1.581887e+09	91.051650	211.127050	2.999696	1.387074	2.848838	3.168573
1.581887e+09	90.844337	211.417287	3.078386	3.123201	3.411111	4.624939
1.581887e+09	90.637025	211.707525	-3.114572	5.431643	4.089597	6.799085
1.581887e+09	90.802875	211.596958	-3.033109	8.097753	4.979262	9.506138
1.581887e+09	91.288579	211.482443	-3.062550	10.840482	6.071373	12.424880
...	...	...	...	...	...	...
1.581888e+09	182.158583	201.452467	-0.986926	0.348276	0.218575	0.411182
1.581888e+09	182.158583	201.397183	-0.978610	0.279135	-0.058413	0.285182
1.581888e+09	182.213867	201.341900	-0.957589	0.193798	-0.283200	0.343162
1.581888e+09	182.158583	201.341900	-0.970083	0.110838	-0.417380	0.431846
1.581888e+09	182.158583	201.286617	-0.936414	0.045190	-0.453966	0.456209

39340 rows × 6 columns

Before linearizing, plotting the head position will help us understand the data.

In [4]:

Copied!





fig, ax = plt.subplots(1, 1, figsize=(10, 10))
ax.plot(
    position_info.head_position_x,
    position_info.head_position_y,
    color="lightgrey",
)
ax.set_xlabel("x-position [cm]", fontsize=18)
ax.set_ylabel("y-position [cm]", fontsize=18)
ax.set_title("Head Position", fontsize=28)
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
ax.plot(
    position_info.head_position_x,
    position_info.head_position_y,
    color="lightgrey",
)
ax.set_xlabel("x-position [cm]", fontsize=18)
ax.set_ylabel("y-position [cm]", fontsize=18)
ax.set_title("Head Position", fontsize=28)

Out[4]:

Text(0.5, 1.0, 'Head Position')

No description has been provided for this image

Specifying the track¶

To linearize, we need a graph of nodes and edges to represent track geometry in the TrackGraph table with four variables:

node_positions (cm): the 2D positions of the graph
edges: how the nodes are connected, as pairs of node indices, labeled by their respective index in node_positions.
linear_edge_order: layout of edges in linear space in order, as tuples.
linear_edge_spacing: spacing between each edge, as either a single number for all gaps or an array with a number for each gap. Gaps may be important for edges not connected in 2D space.

For example, (79.910, 216.720) is the 2D position of node 0 and (183.784, 45.375) is the 2D position of node 8. Edge (0, 8) means there is an edge between node 0 and node 8. Nodes order controls order in 1D space. Edge (0, 1) connects from node 0 to 1. Edge (1, 0) would connect from node 1 to 0, reversing the linear positions for that edge.

For more examples, see this notebook.

In [5]:

Copied!





node_positions = np.array(
    [
        (79.910, 216.720),  # top left well 0
        (132.031, 187.806),  # top middle intersection 1
        (183.718, 217.713),  # top right well 2
        (132.544, 132.158),  # middle intersection 3
        (87.202, 101.397),  # bottom left intersection 4
        (31.340, 126.110),  # middle left well 5
        (180.337, 104.799),  # middle right intersection 6
        (92.693, 42.345),  # bottom left well 7
        (183.784, 45.375),  # bottom right well 8
        (231.338, 136.281),  # middle right well 9
    ]
)

edges = np.array(
    [
        (0, 1),
        (1, 2),
        (1, 3),
        (3, 4),
        (4, 5),
        (3, 6),
        (6, 9),
        (4, 7),
        (6, 8),
    ]
)

linear_edge_order = [
    (3, 6),
    (6, 8),
    (6, 9),
    (3, 1),
    (1, 2),
    (1, 0),
    (3, 4),
    (4, 5),
    (4, 7),
]
linear_edge_spacing = 15
node_positions = np.array(
    [
        (79.910, 216.720),  # top left well 0
        (132.031, 187.806),  # top middle intersection 1
        (183.718, 217.713),  # top right well 2
        (132.544, 132.158),  # middle intersection 3
        (87.202, 101.397),  # bottom left intersection 4
        (31.340, 126.110),  # middle left well 5
        (180.337, 104.799),  # middle right intersection 6
        (92.693, 42.345),  # bottom left well 7
        (183.784, 45.375),  # bottom right well 8
        (231.338, 136.281),  # middle right well 9
    ]
)

edges = np.array(
    [
        (0, 1),
        (1, 2),
        (1, 3),
        (3, 4),
        (4, 5),
        (3, 6),
        (6, 9),
        (4, 7),
        (6, 8),
    ]
)

linear_edge_order = [
    (3, 6),
    (6, 8),
    (6, 9),
    (3, 1),
    (1, 2),
    (1, 0),
    (3, 4),
    (4, 5),
    (4, 7),
]
linear_edge_spacing = 15

With these variables, we then add a track_graph_name and the corresponding environment.

In [6]:

Copied!





sgc.common_position.TrackGraph.insert1(
    {
        "track_graph_name": "6 arm",
        "environment": "6 arm",
        "node_positions": node_positions,
        "edges": edges,
        "linear_edge_order": linear_edge_order,
        "linear_edge_spacing": linear_edge_spacing,
    },
    skip_duplicates=True,
)

graph = sgc.common_position.TrackGraph & {"track_graph_name": "6 arm"}
graph
sgc.common_position.TrackGraph.insert1(
    {
        "track_graph_name": "6 arm",
        "environment": "6 arm",
        "node_positions": node_positions,
        "edges": edges,
        "linear_edge_order": linear_edge_order,
        "linear_edge_spacing": linear_edge_spacing,
    },
    skip_duplicates=True,
)

graph = sgc.common_position.TrackGraph & {"track_graph_name": "6 arm"}
graph

Out[6]:

track_graph_name	environment Type of Environment	node_positions 2D position of track_graph nodes, shape (n_nodes, 2)	edges shape (n_edges, 2)	linear_edge_order order of track graph edges in the linear space, shape (n_edges, 2)	linear_edge_spacing amount of space between edges in the linear space, shape (n_edges,)	linear_edge_specialty denote what edges (denote by 5) are going to be lumped to what edge (denote by 1), shape (n_edges,)
6 arm	6 arm	=BLOB=	=BLOB=	=BLOB=	=BLOB=	=BLOB=

Total: 1

TrackGraph has several methods for visualizing in 2D and 1D space. plot_track_graph plots in 2D to make sure our layout makes sense.

In [7]:

Copied!





fig, ax = plt.subplots(1, 1, figsize=(10, 10))
ax.plot(
    position_info.head_position_x,
    position_info.head_position_y,
    color="lightgrey",
    alpha=0.7,
    zorder=0,
)
ax.set_xlabel("x-position [cm]", fontsize=18)
ax.set_ylabel("y-position [cm]", fontsize=18)
graph.plot_track_graph(ax=ax)
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
ax.plot(
    position_info.head_position_x,
    position_info.head_position_y,
    color="lightgrey",
    alpha=0.7,
    zorder=0,
)
ax.set_xlabel("x-position [cm]", fontsize=18)
ax.set_ylabel("y-position [cm]", fontsize=18)
graph.plot_track_graph(ax=ax)

plot_track_graph_as_1D shows what this looks like in 1D.

In [8]:

Copied!

fig, ax = plt.subplots(1, 1, figsize=(20, 1))
graph.plot_track_graph_as_1D(ax=ax)
fig, ax = plt.subplots(1, 1, figsize=(20, 1))
graph.plot_track_graph_as_1D(ax=ax)

Parameters¶

By default, linearization assigns each 2D position to its nearest point on the track graph. This is then translated into 1D space.

If use_hmm is set to true, a Hidden Markov model is used to assign points. The HMM takes into account the prior position and edge, and can keep the position from suddenly jumping to another. Position jumping like this may occur at intersections or the head position swings closer to a non-target reward well while on another edge.

In [9]:

Copied!





sgc.common_position.LinearizationParameters.insert1(
    {"linearization_param_name": "default"}, skip_duplicates=True
)
sgc.common_position.LinearizationParameters()
sgc.common_position.LinearizationParameters.insert1(
    {"linearization_param_name": "default"}, skip_duplicates=True
)
sgc.common_position.LinearizationParameters()

Out[9]:

linearization_param_name name for this set of parameters	use_hmm use HMM to determine linearization	route_euclidean_distance_scaling How much to prefer route distances between successive time points that are closer to the euclidean distance. Smaller numbers mean the route distance is more likely to be close to the euclidean distance.	sensor_std_dev Uncertainty of position sensor (in cm).	diagonal_bias Biases the transition matrix to prefer the current track segment.
default	0	1.0	5.0	0.5

Total: 1

Linearization¶

With linearization parameters, we specify the position interval we wish to linearize.

In [10]:

Copied!





sgc.common_position.IntervalLinearizationSelection.insert1(
    {
        "position_info_param_name": "default",
        "nwb_file_name": nwb_copy_file_name,
        "interval_list_name": "pos 1 valid times",
        "track_graph_name": "6 arm",
        "linearization_param_name": "default",
    },
    skip_duplicates=True,
)

sgc.common_position.IntervalLinearizationSelection()
sgc.common_position.IntervalLinearizationSelection.insert1(
    {
        "position_info_param_name": "default",
        "nwb_file_name": nwb_copy_file_name,
        "interval_list_name": "pos 1 valid times",
        "track_graph_name": "6 arm",
        "linearization_param_name": "default",
    },
    skip_duplicates=True,
)

sgc.common_position.IntervalLinearizationSelection()

Out[10]:

position_info_param_name name for this set of parameters	nwb_file_name name of the NWB file	interval_list_name descriptive name of this interval list	track_graph_name	linearization_param_name name for this set of parameters
default	chimi20200216_new_.nwb	pos 1 valid times	6 arm	default

Total: 1

And then run linearization by populating IntervalLinearizedPosition.

In [11]:

Copied!

sgc.common_position.IntervalLinearizedPosition().populate()
sgc.common_position.IntervalLinearizedPosition()
sgc.common_position.IntervalLinearizedPosition().populate()
sgc.common_position.IntervalLinearizedPosition()

Computing linear position for: {'position_info_param_name': 'default', 'nwb_file_name': 'chimi20200216_new_.nwb', 'interval_list_name': 'pos 1 valid times', 'track_graph_name': '6 arm', 'linearization_param_name': 'default'}
Writing new NWB file chimi20200216_new_GBGCXYMIWB.nwb

Out[11]:

position_info_param_name name for this set of parameters	nwb_file_name name of the NWB file	interval_list_name descriptive name of this interval list	track_graph_name	linearization_param_name name for this set of parameters	analysis_file_name name of the file	linearized_position_object_id
default	chimi20200216_new_.nwb	pos 1 valid times	6 arm	default	chimi20200216_new_GBGCXYMIWB.nwb	8d132da2-c1e4-402f-ba3a-4b8725a6c87a

Total: 1

Examine data¶

Running fetch1_dataframe will retrieve the linear position data, including...

time: dataframe index
linear_position: 1D linearized position
track_segment_id: index number of the edges given to track graph
projected_{x,y}_position: 2D position projected to the track graph

In [12]:

Copied!





linear_position_df = (
    IntervalLinearizedPosition()
    & {
        "position_info_param_name": "default",
        "nwb_file_name": nwb_copy_file_name,
        "interval_list_name": "pos 1 valid times",
        "track_graph_name": "6 arm",
        "linearization_param_name": "default",
    }
).fetch1_dataframe()
linear_position_df
linear_position_df = (
    IntervalLinearizedPosition()
    & {
        "position_info_param_name": "default",
        "nwb_file_name": nwb_copy_file_name,
        "interval_list_name": "pos 1 valid times",
        "track_graph_name": "6 arm",
        "linearization_param_name": "default",
    }
).fetch1_dataframe()
linear_position_df

Out[12]:

	linear_position	track_segment_id	projected_x_position	projected_y_position
time
1.581887e+09	412.042773	0	90.802281	210.677533
1.581887e+09	412.364853	0	90.520636	210.833775
1.581887e+09	412.686934	0	90.238990	210.990018
1.581887e+09	412.488270	0	90.412714	210.893645
1.581887e+09	412.007991	0	90.832697	210.660660
...	...	...	...	...
1.581888e+09	340.401589	1	175.500994	212.958497
1.581888e+09	340.373902	1	175.477029	212.944630
1.581888e+09	340.394065	1	175.494481	212.954729
1.581888e+09	340.346214	1	175.453064	212.930764
1.581888e+09	340.318527	1	175.429100	212.916898

39340 rows × 4 columns

We'll plot the 1D position over time, colored by edge, and use the 1D track graph layout on the y-axis.

In [13]:

Copied!





fig, ax = plt.subplots(figsize=(20, 13))
ax.scatter(
    linear_position_df.index,
    linear_position_df.linear_position,
    c=linear_position_df.track_segment_id,
    s=1,
)
graph.plot_track_graph_as_1D(
    ax=ax, axis="y", other_axis_start=linear_position_df.index[-1] + 10
)

ax.set_xlabel("Time [s]", fontsize=18)
ax.set_ylabel("Linear Position [cm]", fontsize=18)
ax.set_title("Linear Position", fontsize=28)
fig, ax = plt.subplots(figsize=(20, 13))
ax.scatter(
    linear_position_df.index,
    linear_position_df.linear_position,
    c=linear_position_df.track_segment_id,
    s=1,
)
graph.plot_track_graph_as_1D(
    ax=ax, axis="y", other_axis_start=linear_position_df.index[-1] + 10
)

ax.set_xlabel("Time [s]", fontsize=18)
ax.set_ylabel("Linear Position [cm]", fontsize=18)
ax.set_title("Linear Position", fontsize=28)

Out[13]:

Text(0.5, 1.0, 'Linear Position')

We can also plot the 2D position projected to the track graph

In [14]:

Copied!





fig, ax = plt.subplots(1, 1, figsize=(10, 10))
ax.plot(
    position_info.head_position_x,
    position_info.head_position_y,
    color="lightgrey",
    alpha=0.7,
    zorder=0,
)
ax.set_xlabel("x-position [cm]", fontsize=18)
ax.set_ylabel("y-position [cm]", fontsize=18)
ax.plot(
    linear_position_df.projected_x_position,
    linear_position_df.projected_y_position,
)
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
ax.plot(
    position_info.head_position_x,
    position_info.head_position_y,
    color="lightgrey",
    alpha=0.7,
    zorder=0,
)
ax.set_xlabel("x-position [cm]", fontsize=18)
ax.set_ylabel("y-position [cm]", fontsize=18)
ax.plot(
    linear_position_df.projected_x_position,
    linear_position_df.projected_y_position,
)

Out[14]:

[<matplotlib.lines.Line2D at 0x7f5d0807e760>]

Interactive selection¶

Note: Work in Progress

`NodePicker`¶

Linearization heavily depends on the track graph is specified. To help simplify setting the nodes/edges, we can use the NodePicker to interactively set node positions and edges based on the video.

In [15]:

Copied!





%matplotlib widget

key = {
    "nwb_file_name": nwb_copy_file_name,
    "interval_list_name": "pos 1 valid times",
}

epoch = (
    int(
        key["interval_list_name"]
        .replace("pos ", "")
        .replace(" valid times", "")
    )
    + 1
)
video_info = (
    sgc.common_behav.VideoFile()
    & {"nwb_file_name": key["nwb_file_name"], "epoch": epoch}
).fetch1()

io = pynwb.NWBHDF5IO("/stelmo/nwb/raw/" + video_info["nwb_file_name"], "r")
nwb_file = io.read()
nwb_video = nwb_file.objects[video_info["video_file_object_id"]]
video_filename = nwb_video.external_file.value[0]

fig, ax = plt.subplots(figsize=(8, 8))
picker = sgc.common_position.NodePicker(ax=ax, video_filename=video_filename)
%matplotlib widget

key = {
    "nwb_file_name": nwb_copy_file_name,
    "interval_list_name": "pos 1 valid times",
}

epoch = (
    int(
        key["interval_list_name"]
        .replace("pos ", "")
        .replace(" valid times", "")
    )
    + 1
)
video_info = (
    sgc.common_behav.VideoFile()
    & {"nwb_file_name": key["nwb_file_name"], "epoch": epoch}
).fetch1()

io = pynwb.NWBHDF5IO("/stelmo/nwb/raw/" + video_info["nwb_file_name"], "r")
nwb_file = io.read()
nwb_video = nwb_file.objects[video_info["video_file_object_id"]]
video_filename = nwb_video.external_file.value[0]

fig, ax = plt.subplots(figsize=(8, 8))
picker = sgc.common_position.NodePicker(ax=ax, video_filename=video_filename)

Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …

After selection, we can retrieve the data using the node_positions and edges attributes

In [16]:

Copied!

picker.node_positions
picker.node_positions

Out[16]:

array([], dtype=float64)

In [17]:

Copied!

picker.edges
picker.edges

Out[17]:

[[]]

Selector¶

We can also draw a 2d polygon around the track and attempt to recover the graph.

In [31]:

Copied!





%matplotlib widget

fig, ax = plt.subplots(figsize=(8, 8))
selector = sgc.common_position.SelectFromCollection(ax, video_filename)

print("Select points in the figure by enclosing them within a polygon.")
print("Press the 'esc' key to start a new polygon.")
print("Try holding the 'shift' key to move all of the vertices.")
print("Try holding the 'ctrl' key to move a single vertex.")
%matplotlib widget

fig, ax = plt.subplots(figsize=(8, 8))
selector = sgc.common_position.SelectFromCollection(ax, video_filename)

print("Select points in the figure by enclosing them within a polygon.")
print("Press the 'esc' key to start a new polygon.")
print("Try holding the 'shift' key to move all of the vertices.")
print("Try holding the 'ctrl' key to move a single vertex.")

Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …

Select points in the figure by enclosing them within a polygon.
Press the 'esc' key to start a new polygon.
Try holding the 'shift' key to move all of the vertices.
Try holding the 'ctrl' key to move a single vertex.

Up Next¶

Next, we'll combine ephys and position data in a process called ripple detection.

Position - Linearization¶

Overview¶

Imports¶

Retrieve 2D position¶

Specifying the track¶

Parameters¶

Linearization¶

Examine data¶

Interactive selection¶

NodePicker¶

Selector¶

Up Next¶

`NodePicker`¶